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Ensemble averaged dynamic modeling 

By D. Carati 1 , A. Wray 2 AND W. Cabot 3 


The possibility of using the information from simultaneous equivalent large eddy 
simulations for improving the subgrid scale modeling is investigated. An ensemble 
average dynamic model is proposed as an alternative to the usual spatial average 
versions. It is shown to be suitable independently of the existence of any homogene- 
ity directions, and its formulation is thus universal. The ensemble average dynamic 
model is shown to give very encouraging results for as few as 16 simultaneous LES's. 


1. Introduction 

The equation for large eddy simulation (LES) is obtained by applying a spatial 
filter to the Navier-Stokes equation. The LES equation thus describes the evolution 
of a filtered velocity field Hi which explicitly depends on the small scales through 
the subgrid scale stress r t j = UiUj — Ui uy. 

d t Ui +djUjUi = -djT t j. ( 1 . 1 ) 

For simplicity, we only consider incompressible flows. The pressure p is then chosen 
to satisfy the incompressibility condition. Clearly, is a large scale quantity 
depending mainly on the small scale velocity field. However, it is usually modeled 
as a function of the resolved velocity field as in the Smagorinsky eddy viscosity 
model (Smagorinsky, 1963): 

r i} - ^T kk 6ij * -2CA 2 |5|5 (> , (1.2) 

where S,y = \ (d t Uj + djUi) and \S\ = (2 S x jS x j) l l 2 . In the original formulation 
of the Smagorinsky model, the parameter C must be obtained from some fitting 
procedure. Recently, this model has been improved by the introduction of the 
dynamic procedure, which allows a self calibration of the parameter C and gives 
an explicit expression as a function of the resolved field C = C(ufc)- However, any 
procedure that determines the subgrid scale stress in terms of the resolved field can 
only be an approximation. Indeed, the same resolved field may be compatible with 
many different small scale velocity fields. This is reflected in the a priori tests which 
show very poor correlations between the models r,j » r t ^(u*) and the actual r tJ 
obtained from direct numerical simulations (see Winckelmans et a/, in this volume). 
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Clearly more information is needed to properly reconstruct the subgrid-scale 
stress. The introduction of stochastic model for r tJ is a first attempt to intro- 
duce models that are not fully determined by the resolved field (Carati et al , 1995; 
Chasnov, 1991; Leith, 1990; Mason & Thomson, 1992). Here, we explore another 
approach which consists in running simultaneously several statistically equivalent 
LES’s and constructing the model by using information from the set of resolved 
velocity fields: 

d t u T { + djU’-tf = -d t p r + vo v 2 < -djTlj r = 1, . . . , R. ( 1.3) A 

Here, r is a new index corresponding to the realization and R is the total number of 
realizations. The concept of statistically equivalent LES’s will be defined in Section 
3. The model we propose to test should generalize the classical subgrid scale model 
(r£ = r[j(u r k )) by allowing an explicit dependence on the velocity field from other 
members in the set: 


r[ } = r^aui}) ( 1 . 4 ) 

Clearly, in that case the subgrid scale model in the LES labeled r will not be a 
function of the resolved velocity field u r k only. 

In the following section, we will present the dynamic procedure and its general- 
ization to several LES’s. We also present an alternative formalism to the classical 
dynamic model. Some results for decaying and forced isotropic turbulence and for 
channel flow are discussed. 

2. The dynamic procedure 

The dynamic procedure is based on an exact relation between subgrid scale 
stresses for different filter widths (Germano, 1992; Ghosal et al 1995; Lilly, 1992). 
This relation is obtained by introducing a second filter Gt , usually referred to as 
the test filter , denoted by we will call the original filter G\. The application of 
this new filter to Eq. (1.1) yields: 

d t u, + djUjUi = -dip + voV 2 u, - djT i} - d]L, } , ( 2 . 1 ) 

where = u x \ij — u t u ; is the Leonard tensor. This equation governs the evolution 
of the field u obtained by the application of the filter G 2 = G t ★ G 1 to the fully 
resolved velocity. Thus, an equivalent equation should be obtained by applying G 2 
directly to the Navier-Stokes equation: 

d t Uj + djUjUi = -dip + i'o V 2 U, - djTij. (2.2) 

Here, the subgrid stress tensor is defined by T, ; = Ui~uj — u, Uj. The comparison 
between equations (2.1) and (2.2) readily leads to the Germano identity: 


Lij + ?ij - Tij = 0 . 


(2.3) 
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When approximate models r tJ « rfj and Tij as are used, this identity is 

violated. However, the error E tJ = Lij + t™ — Tff ^ 0 may be used to calibrate 
the models. When the Smagorinsky model is used at both grid and test levels, the 
error is a linear function of the Smagorinsky parameter: 

Eij = + Cfiij - Cotij (2.4) 


where 


Qij = -2A \S\Stj 
Pa = -2A 2 \S\S tj 

The calibration of C is usually achieved by using a least square method for mini- 
mizing Eij. The integral 


I[C} = 




(y) 


(2.5) 


is thus minimized with respect to C . 

A first difficulty encountered when using the dynamic procedure for determining 
C has been pointed out by Ghosal et al (1993,1995), who showed that this procedure 
requires the solution of an integral equation for C unless both of the following 
conditions are satisfied: 

1. There are one or more directions of homogeneity in the flow. 

2. The flow is fully resolved in the other directions ). 

In that case, C is assumed to be constant along the direction of homogeneity and 
can be taken out of the test filter operation'"'. Moreover, the flow being fully resolved 
in the other direction(s), the test filter must only act in the homogeneous direction. 
The error (2.4) then reduces to: 


Eij — Lij + CMij ( 2 . 6 ) 

where M tJ = a l} — (j tJ and the dynamic prediction for C reads: 

q _ {LjjMjj)h 
” (MijMij)h 

where the brackets ( )/, represent a spatial average in the homogeneous direction(s). 
If the two aforementioned conditions for replacing expression (2.4) by (2.6) are not 
fulfilled, one could argue that C is slowly varying in space and that (2.6) should be 
a valid approximation independently of the existence of a direction of homogeneity. 
The minimization of the global quantity I[C] then leads to a local expression for C: 
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Unfortunately, this approximation has proved to he very poor, and the result- 
ing C depends strongly on space. Since in almost all LES's at least one of the 
aforementioned conditions is violated, a mathematically clean implementation of 
the dynamic model always requires the solution of an integral equation (Ghosal et 
al 1995) . 

A second difficulty with the dynamic model is that C takes negative as well as 
positive values. Positive values correspond to the classical eddy dissipation picture 
for the subgrid scales. The negative values were first interpreted as the capability 
of the dynamic model to predict reverse energy transfer (backseat ter). Unfortu- 
nately, the modeling of backseat ter by a negative Smagorinsky coefficient leads to 
numerical instabilities. This problem is easily solved by constraining a priori the 
minimization of I[C] so that only positive values of C are accepted. The resulting 
C (obtained either by solving an integral equation or by using a spatial average) is 
the same as before but clipped to positive value. Thus, C must then be replaced 
by (C + \C\)/2. Although this clipping procedure can be derived properly from a 
constrained minimization procedure, it is usually considered an undesirable exten- 
sion of the dynamic model. In particular, the clipping corresponds to turning off 
the model when' the dynamic procedure "tries to build a model for backseat ter.” 
In some sense, the resulting model does not use all the information available from 
the dynamic procedure. Hence, it is desirable to have a dynamic model with as few 
clipped values as possible for C. 

We will discuss in the following sections how the simultaneous use of several 
statistically equivalent LES’s may solve these two difficulties. 

3 Statistical LES & dynamic model 

3.1 Definition of the ensemble 

We first, discuss the problem of defining the ensemble of runs needed for the 
statistical tests without considering the modeling problem. The equations (1.3) 
correspond to H different LES’s. In order to have a "good” ensemble, these LES’s 
should correspond to statistically equivalent and statistically independent realiza- 
tions of the same problem. Although these requirements art? intuitively clear, it is 
worthwhile to define them as properly as possible. The first step consists in defining 
precisely what is an "acceptable” simulation for a given problem. From the strict 
mathematical point of view, a flow described by the Navier-Stokes equation or by 
an LES equation is completely defined by the knowledge of 

1. The domain T> in which the flow is considered. 

2. The conditions on the boundary dV of this domain v(dT>,t) = /(f). 

3. The initial conditions e(x, 0) = cq ( x ) Vr £ V. 

However, in a simulation of a turbulent flow only the domain and the boundary 
conditions are rigorously fixed. Indeed, because of the 1 lack of sensitivity to initial 
conditions in a turbulent flow, different simulations with different initial conditions 
sharing some properties are considered to characterize the same flow. Thus, the 
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requirement that the initial conditions are known is usually replaced by some weaker 
constraints, and point (3) is replaced by 

3’. The initial conditions t>(x,0) = v 0 (x; wi) are generated using random num- 
bers wi and satisfy some constraints: P^o] = s = 1, ... S'. 

For example, in homogeneous turbulence, the first constraint s = 1 will be on 
the spectrum of Vq. For the channel flow, one could impose the profile of the 
velocity and of the fluctuation in each direction. We will not discuss in detail the 
minimal constraints that must be imposed on the initial conditions in order to have 
a reasonable simulation. We only suppose that these constraints do exist. Now, it 
is possible to give some precise definition of the ensemble of LES’s: 

Definition 1: Two LES's are statistically equivalent if the domain of the flow 
and the boundary conditions are exactly the same and if the initial 
conditions satisfy the same set of constraints. 

Definition 2: Two LES's are statistically independent if the initial conditions 
are generated with uncorrelated random numbers w\. 

For a stationary flow, such equivalent and independent initial conditions can be 
obtained by running a single LES and recording several velocity fields separated by 
at least one large eddy turnover time when turbulence is fully developed. 

3.2 Ensemble average dynamic model 

In what follows, we will focus on a simple generalization of the Smagorinsky 
model which reads: 


* -2CA 2 |S r |S: r (3.2.1) 

Thus, we basically use the Smagorinsky model in every realization with the following 
additional assumption : 

Hypothesis 1 : The Smagorinsky coefficient is independent of the realization for sta- 
tistically equivalent flows. 

This assumption defines the model in such a way that the unknown parameter in 
the LES is “universal”. The formulation thus mixes some aspects of both LES and 
Reynolds average simulations. 

The dynamic procedure can also be used to determine C when several LES’s are run 
in parallel. In that case, the model depends on the resolved flow from other real- 
izations (1.3). Indeed, the quantity that needs to be minimized is a straightforward 
generalization of I[C]: 

AC] = £ I <*y£(^(y)) 2 < 3 - 2 - 2 > 

r J V if 

where now E x j as well as L XJ , j3 tJ , and a t j depend on the realization (EJj ~ L \ + 
C j3fj — Ca^). We now make another assumption: 
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Hypothesis 2: For large ensembles , the Smagorinsky coefficient is slowly dependent 
on space and can be taken out of the test filter. 

The quantity 2[C] then reduces to 

r= 1 i j 

which leads to the same expression for C as in the spatial average version of the 
dynamic procedure: 


r = (LvMnl 

where now the brackets represent an ensemble average. We will see in the next 
section that hypothesis 2 is very well justified by the numerical results. 

3.3 Alternative formalism for the dynamic model 

The usual formulation of LES Eq. (1.1) is not fully satisfactory because the evo- 
lution of the filtered velocity is given in terms of quantities that are not filtered, 
whereas all numerically computed quantities are filtered in some way. This is well 
known, but, to our knowledge, its effect on the dynamic model formulation has 
never been carefully considered. In this section, we propose an alternative dynamic 
model formulation which should be fully self-consistent with the filtered equation 
for the resolved field. First, we assume that all the quantities in the LES equation 
are filtered and Eq. (1.1) must then be replaced by 

d t u t — vq V 2 TT, — djtljUi — djT tJ — d t p. (3.3.1) 

This redefines the subgrid scale stress as 


T tJ = UiUj ~ U t Uj 

The application of the test filter to the LES equation (3.3.1) yields: 


d t u, + dj(ujUi) = v 0 V 2 Ui - d t p - djT i} - djL,j , 

and the comparison with the “one-step” application of Gz to u t leads to the following 
equality: 

Lij + t ij - Tij = 0. (3.3.2) 


where now 
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At this point it is important to ensure that the model for the subgrid scale is 
also expressed in terms of a filtered quantity. The simplest generalization of the 

Smagorinsky model would then be rj] — Cf3 t j and T tJ — Ca x) . The dynamic 
procedure is then easily implemented and yields 


c _ (LijNji) 

(Fifia) 


(3.3.3) 


where N tJ = fi tJ — Of course, the expression (3.3.3) also relies on the assumption 
that C can be taken out from the filtering operators. This assumption is very 
important here because, in the equality (3.3.2), the Smagorinsky coefficient only 
appears in filtered quantities. This means that the integral equation formulation of 
this alternative dynamic model would be much more complicated than the classical 
formalism. However, if hypothesis 2 is valid, the present formalism appears to be 
more consistent with the LES equation. 


4. Test on isotropic turbulence 

4-1 Decaying turbulence 

The statistical average dynamic model described in section 3.1 has been tested 
in decaying turbulence for 32 3 LES. A first series of numerical experiments have 
determined how large the ensemble of simultaneous LES's must be ( i.e . how large 
R should be). The criteria used to determine the minimal size of the ensemble were 
focused on 

1. The spatial variability of C. 

2. The percentage of negative C. 

3. Comparison with the volume average dynamic model. 

4. Comparison with direct numerical simulations. 

The first conclusion we have reached is quite encouraging. Indeed, it appears that 
with only 16 simultaneous LES’s, the ensemble average dynamic model performs as 
well as the volume average model. The spatial variability of C decreases drastically 
when R increases (see Fig. 1). This is also reflected on the probability distribution 
function (PDF) of C (see Fig. 2). 

The comparison between a 512 3 DNS and dynamic model shows good agreement 
both for the total resolved energy (see Fig. 3) and for the spectra. The results 
for R = 16 are indistinguishable for the volume average and for the ensemble aver- 
age. Here the comparison with the dynamic model has been made by running an 
ensemble of unrelated volume average LES’s. This allows comparison of the both 
the means and the standard deviations. The standard deviations are computed for 
the 3-d energy spectra at each A', and quantities such as total resolved energy and 
compensated spectra are then computed from the mean and mean±a spectra. 
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FIGURE 1. Typical profile of C in decaying isotropic turbulence. R = l: — • — ; 

R=4:— ♦ — ; R=16: . 



FIGURE 2. PDF of C in decaying isotropic turbulence. Symbols as in Fig. 1. 

4.2 Forced turbulence 

We have run an ensemble of 32 3 forced turbulence LES's with zero molecular 
viscosity. Fig. 4 shows that the mean resolved energy and the standard deviation 
evolve in a very similar way for both the volume and the ensemble average models. 
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FIGURE 3. Energy decay: comparison with DNS and volume average. DNS: ; 

ensemble- averaged (mean):n ; ensemble-averaged (mean 4-sigma): ; ensemble- 

averaged (mean-sigma): ; volume-averaged (mean):* : volume- averaged 

(mean 4-sigma): ; volume- averaged (mean-sigma): 



FIGURE 4. Resolved energy in forced isotropic turbulence: average vs volume. 
Symbols as in Fig. 3, without DNS. 




246 


D . Carati 7 A. Wray & W. Cabot 


4.0 


3.0 


0 ) 

c 

a> 

r> 

Q) 

« 

00 

c 

a) 

a. 

E 

o 

O 


2.0 


1.0 


/ 


So,'/ 

M // - 


n 
; / 


/ / 

4 


o.o 


0.0 


5.0 


10.0 


15.0 


FIGURE 5. Compensated energy spectrum in forced isotropic turbulence: average 
vs volume. Symbols as in Fig. 3. 

It is also very important to notice that the standard deviation saturates so that the 
LES's in the ensemble do not evolve towards very different states. 

It is also interesting to compare the compensated energy spectrum to check if an 
inertial range is observed. Of course, with 32* LES, we do not expect to obtain 
a very good estimate of the Kolmogorov constant. However, the results in Fig. 5 
show that the observed "Kolmogorov constant” is in a reasonable range of values. 
These spectra are at time ~ 27 in the units of Fig. 4. 

5. Tests in channel flow 

We did not reach the stage of “production runs” for the channel flow, so the 
tests presented here are very preliminary and have been focused on the behavior 
of C as a function of the ensemble size (i?). The data collected from the runs 
concern the PDF of C and the fraction of negative C . Because of the channel flow 
inhomogeneity, the PDF of C depends on the wall normal coordinate. However, 
the trends for increasing numbers of runs (R) is similar across the channel, and we 
only present in Fig. 6 the results for y = 0.1 near midchannel. 

We also show the fraction of negative C (Fig. 7). Since the channel flow simula- 
tions used in these tests have a non-zero molecular viscosity, the relevant stability 
condition is the percentage of C leading to a total (molecular + eddy) negative 
viscosity. Here again, the results are encouraging for /? ~ 1C (less than 15% of the 
points need to be clipped) while the local version of the dynamic model for only 
one LES requires about 40% clipping. 
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FIGURE 6. Probability distribution function of C for different ensemble sizes at 
y = 0.1. R=l:* ; R=2:« ; R=4: + ; R=8;o ; R=16:n . 



FIGURE 7. Fraction of negative total viscosity as a function of y for different en- 
semble sizes. Symbols as in Fig. 6. 


6. Conclusion 

The statistical tests presented in this report have shown that the knowledge of 
statistically equivalent resolved velocity fields may be useful in deriving new subgrid 
scale models. We have used the additional information available from the different 
LES’s to create an ensemble average version of the dynamic model. This dynamic 
model has the following advantages: 
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1. A local version of the ensemble average dynamic model may be derived in the 
limit of large ensemble sets. 

2. The local formulation does not rely on any homogeneity assumption. It can thus 
be adapted to any geometry, unlike to the classical volume (or surface or line) 
average dynamic model. 

3. The theoretical limit of large ensemble sets is closely approached for R « 16. 
This is indicated by the PDF of C, which is very peaked for R — 16. Also, the 
spatial variations of C decrease drastically for increasing ensemble sizes and seem 
to be quite mild for R = 16. 

For the examples treated in this work (decaying and forced isotropic turbulence), 
the volume average version of the dynamic model is justified. Remarkably, in those 
cases, the results from the ensemble average and the volume average versions are 
almost indistinguishable. 

The next interesting step in the investigation of statistical LES is to apply this 
model to fully inhomogeneous problems (for which the mathematically consistent 
classical dynamic model requires the solution of an integral equation). The addi- 
tional cost of multiple simultaneous LES’s may be ameliorated by a reduction in 
the time of simulation since the statistics should converge more rapidly. 
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